Saturation mechanism of the Weibel instability in weakly 

magnetized plasmas 

Tsunehiko N. KatcQ 

National Astronomical Observatory of Japan, 
2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan 
(Dated: 27 June, 2005) 

Abstract 

The saturation mechanism of the Weibel instability is investigated theoretically by considering 
the evolution of currents in numerous cylindrical beams that are generated in the initial stage of 
the instability. Based on a physical model of the beams, it is shown that the magnetic field strength 
attains a maximum value when the currents in the beams evolve into the Alfven current and that 
there exist two saturation regimes. The theoretical prediction of the magnetic field strength at 
saturation is in good agreement with the results of two-dimensional particle-in-cell simulations for 
a wide range of initial anisotropy. 



*Electronic address: katoutn@cc.nao.ac.jp 



1 



Collisionless plasmas with anisotropic velocity distributions drive the Weibel instability 
, thereby generating magnetic fields. Recently, this instability in astrophysical as well 
as laboratory plasmas has attracted considerable attention. For example, it is considered 
that this instability can be driven in strong collisionless shock waves associated with various 
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astrophysical phenomena, e.g., pulsar winds [3|,|4|, gamma-ray bursts and/or their afterglows 
0| , or gravitational collapse of large-scale structures in the universe It can also be driven 
if temperature gradients are present in plasmas [7]. The magnetic field generated by the 
instability is responsible for synchrotron and/or "jitter" radiation from the existing high- 
energy particles jq]. Furthermore, such magnetic fields can provide an effective scattering 
mechanism for charged particles. For example, it would affect the dissipation process of 
collisionless shock waves, efficiency of the Fermi acceleration, or heat conductivity by charged 
particles. In all these cases, the amplitude of the magnetic field is of primary importance. 

The magnetic field strength attains the maximum va 
Several conditions for saturation have been proposed 
may be unsatisfactory because the role of current was not treated explicitly in these studies. 
Since the instability enters the nonlinear regime before saturation, the currents and fields 
must be considered self-consistently. As shown by numerical simulations, the currents are 
carried by numerous cylindrical beams 11, 12, Therefore, the physics of saturation will 



ue at the saturation of the instability. 
3 0, [l^ . However, these conditions 



be elucidated by examining the evolution of such beams. 

In this letter, the nonlinear saturation mechanism of the Weibel instability, i.e., that of 
the transverse modes of filamentation instability, is investigated with a physical model of the 
beams. The theoretical analysis and simulations are performed within a two-dimensional 
framework. (Therefore, in three dimensions, the results might be modified by other effects.) 
The method is applicable to electron-proton as well as electron-positron plasmas. However, 
for electron-proton plasmas, only the saturation due to the electron currents is considered. 
Although proton currents later cause the second saturation [3], it is beyond the scope of the 
present study and will be investigated in the future. In the following analysis, protons are 
treated clS cL charge neutralizing background. 

Theoretical model of saturation — In this section, the saturation mechanism is modeled 
by using Gaussian units and the following notation: speed of light, c; velocity of a parti- 
cle normalized by the speed of light, j3; Lorentz factor of the particle, 7 = (1 — /3 2 ) -1 / 2 ; 
electron mass, m e ; electron charge, — e; proton mass, m v \ mean electron number density, 
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n e ; electron plasma frequency, u) pe = (47rn e e 2 /m e ) 1 / 2 ; electron skin depth, l = c/u pe « 
5.3 x 10 5 n~ 1 / 2 [cm]. 

The initial condition is set as follows. The plasma is initially uniform and unmagnetized 
(or weakly magnetized) with an anisotropic velocity distribution. The initial velocity dis- 
tributions of the particles are axisymmetric, and the velocity dispersion along the z-axis 
is larger than those along the other axes. For electron-positron plasmas, the distribution 
functions are the same for both species. For electron-proton plasmas, the distribution func- 
tion of only the electrons is considered. Such a distribution function would be plausible for 
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two- stream like conditions 



13, 



The basic physical mechanism of the instability is described as follows. In the initial 
stage of the instability, the charged particles in the plasma are deflected by small magnetic 
fluctuations. In the present condition, most particles have larger velocities along the z- 
axis. Since the particles carrying current toward the +z direction and — z direction are 
deflected in opposite directions, they are separated into different regions, thereby producing 
net currents in the plasma that generate magnetic fields 5] . As the amplification of the 
magnetic fields increases, the distance between the two populations of particles increases, 
and the currents are amplified, and vice-versa. Thus, both currents and magnetic fields 
increase exponentially. In this evolution, it is important for the instability to enter the 



non 
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inear regime before saturation. As revealed by two- or three-dimensional simulations 
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13j , many isolated cylindrical beams (or current filaments) are formed in the plasma 
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121 ] even before saturation. 



long before saturation; each beam carries a net current and generates a magnetic field around 
itself. Since two beams with currents in the same direction attract each other and tend to 
coalesce into a larger beam, the beams grow with time 

Based on these observations, we model the nonlinear evolution of the instability from the 
"initial condition" in which the plasma consists of many cylindrical beams to saturation. 
Initially, the radius and magnitude of current in the beams are almost uniform because the 
plasma is homogeneous. Each beam carries a net current in either the +z or —z direction. 
Subsequently, both the current and magnetic field in each beam increase exponentially due 
to particle separation in the beam. The radii of the beams also increase due to the coalescing 
process, although it is considerably slower than the increase in current and magnetic field. 
However, they are still considered to be uniform because they essentially evolve equally. 
Therefore, the plasma can be modeled as an ensemble of uniform cylindrical beams with 
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the same radius r and magnitude of current /. (A similar model was recently considered in 
for investigating the evolution after saturation.) By assuming a uniform current 
density J, where I = irr 2 J, the root-mean-square value of the magnetic field strength within 
a cylindrical beam is calculated using the Ampere's law as follows: 

B = V2nrJ/c = V2I/(rc). (1) 

Although in practice, the magnetic field is affected by neighboring beams, their net effect 
will be low and the order of the magnetic field strength in the beam will be given by Eq. (JJJ) . 
Thus, we assume that the magnetic field (JIJ in each beam is essential to the saturation 
process and those generated by the neighboring beams are important only in the coalescing 
process. 

On the basis of this model, we consider the evolution of one of the beams. The exponential 
increase in current and magnetic field ceases when the magnitude of the current is equal to 

n 

the smaller of the following two upper limits. One limit is the Alfven current [la], which 
occurs due to the self-generated magnetic field; it is expressed in Gaussian units by 

lA = hW, (2) 

where Jo = m e c 3 /e (~ 17000 [A]), /3n is the magnitude of the z component of f3, and the 
angle brackets ( ) denote an average taken over the beam volume. (This limit was originally 
derived for a monoenergetic single-directed particle beam. In this case, the averaged value 
of 7/3|| is used instead of that in the monoenergetic cases because the energy of particles 
is distributed.) We can apply this limit even when the separation process is in progress 
and a fraction of particles carry current in the opposite direction because the self-generated 
magnetic field is associated with the net current I, and it always ensures that I does not 
exceed I a- The other limit is the maximum current Ip that can be carried by all the particles 
within the beam: 

Ip = nr 2 Jp, where Jp = efxn e (/3\\)c (3) 

with [l = 1 for electron-proton plasmas and \x = 2 for electron-positron plasmas. This limit 
is reached when the separation of the current density is completed before I equals I a- Since 
the distribution of particles in the beam approaches an isotropic distribution depending on 
the evolution of the instability, we must consider this effect; it is expressed by a factor \ as 
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follows: 

((3 l{ ) = X {/3\\)o, W = xWo, (4) 

where subscript "0" denotes an initial value. Since a complete isotropization of the strong 
initial anisotropy yields the lower limit x = 1/3, we observe that 1/3 < x < 1- 

After the exponential growth ceases, the beam evolves through the coalescing process. 
In the former case (I a < Ip), which is termed the Alfven limit, once the current equals the 
Alfven current I a, it no longer increases and it retains the value of I a, while the radius 
continues to increase. Therefore, from Eq. (fl)l. the magnetic field decreases monotonically 
thereafter. In the latter case (I a > Ip), which is termed the particle limit, current gradually 
increases until it reaches the Alfven current; during this period, it increases with the radius 
while satisfying the condition / = nr 2 Jp, and the magnetic field also increases. After the 
current equals the Alfven current, the magnetic field decreases monotonically, as in the 
Alfven limit. 

Therefore, for both limits, the magnetic field strength becomes maximum when the cur- 
rent evolves into the Alfven current. From Eq. the maximum value -B max is expressed 
as follows: 

B max = V2I A /(fc) = v/^uHf/Jo)- 1 ^, (5) 

where f is the beam radius at saturation and B* is the magnetic field strength given by 
= c(Aim e m e ) l l 2 m 3.2 x 10~ 3 rig' 2 [G]. In the Alfven limit, the radius at saturation ta 
is related to the initial radius vq or the wavelength of the most unstable mode in the linear 
theory Ao- Here, it is simply expressed using two factors a and ( as follows: 

f A = ar = a(\ . (6) 

Thus, the magnetic field at saturation is given by 

B max = V2c a (jP\\}o (Xo/ky 1 £*, (7) 

where ca = XA/{ a C)i an d Xa denotes the effect of isotropization at saturation [see Eq. (j3J)]. 
In the particle limit, the radius at saturation fp is determined from the condition nf P Jp = Ia 
as follows: 

f P = (I A /irJ P ) 1/2 = 2/ ((7/5||)o//i(A)o) 1/2 (8) 
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which is on the order of Iq for nonrelativistic cases. Using this radius, the maximum magnetic 
field strength is given by 

Smax = XP (/i(7/?||)o(A|) /2) 1/2 £*, (9) 

where xp is the isotropization factor at saturation. As mentioned earlier, the values of 
-Bmax in Eqs. (jjj) and can be regarded as typical magnetic fields of the entire plasma at 
saturation in the respective cases. The two parameters (ca and xp) i n these expressions are 
expected to be of order unity. 

The ratio of the maximum magnetic energy density to the initial particle kinetic energy 
density, 77, is written as 77 = (B max / B*) 2 [2[/,((y) — (It should be noted that for 

electron-proton plasmas, the kinetic energy of protons is not considered in this equation.) 
For the particle limit regime, 77 yields a subequipartition value. In particular, rj ~ Xp/2 for 
nonrelativistic plasmas and 77 ~ Xp/4 for ultrarelativistic plasmas. The conclusion that 77 
yields a subequipartition value for strongly anisotropic cases is in agreement with that of 



previous simulations 



Comparison with numerical simulations — The results of the developed analytical model 
are compared with those of numerical simulations for electron-positron and electron-proton 
plasmas (m p /m e = 1836). The simulation code used is a relativistic, electromagnetic, 
particle-in-cell code with two spatial and three velocity dimensions. This code was de- 
veloped based on a general description by Ref. [17]. The x-y plane perpendicular to the 
z-axis is considered to be the simulation plane. The initial particle distribution is expressed 
in terms of the normalized momentum u = 7/3 that is common for all the species. Each 
component of u obeys the Gaussian distribution with a standard deviation of cry for the 
z component or a± for the other components. The simulations were performed using a 
512 x 512 grid with ~ 50 particles per cell per species under periodic boundary conditions 
for several values of cry ranging from 0.12 to 10 with fixed aj_ = 0.1. For each simulation, the 
physical size of the simulation box in each direction, L, was considered to be at least 7 times 
larger than the typical beam diameter at saturation, which was estimated by a preliminary 
simulation [e.g., L = 120Zo for an = 0.12 (largest); L = 30/o for cry = 0.6 (smallest)]. Then, 
in each simulation, the evolution of the magnetic field strength averaged over the simulation 
box was monitored and its maximum value was obtained; the calculation time T considered 
was long enough to obtain the maximum value [e.g., T = 1500^,3 for cry = 0.12 (longest); 



T = 25ui~e for a\\ = 1.1 (shortest)]. Figure 1 shows the obtained maximum values denoted 
by dots as a function of the initial anisotropy, <J\\/<7± — 1: (a) electron-positron plasma and 
(b) electron-proton plasma. 

The corresponding theoretical results are obtained as follows: First, in Eqs. (J7J) and 
we approximate that (7/3||)o = cry an d (/3||)o = cr ||/7> where 7 = (1 + cr| + 2a 2 J 1 ^ 2 . Next, 
in the present case, since the linear theory of nonrelativistic weak anisotropic plasmas 
applicable to the Alfven limit regime, we obtain 

-1/2 
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A = 27r^3//i[(a,|/a ± ) 2 -lj l . (10) 

(It should be noted that, when a± > 1 or a\\ > 1, the relativistic dispersion relation must 
be employed.) Finally, using these expressions, we obtain 

^= = {yf £ 1rh/^ 2 - 1 ] 1/2 (Aifv& Hmi t ) (n) 

XP&W [/^/ (2 ; y)] 1 ^ 2 (particle limit). 



The results are plotted in Fig. 1; dotted curves represent the results for the Alfven limit and 
solid curves for the particle limit. The parameters (ca and xp) are assumed to be constants. 
Their values are taken as {ca,Xp) = (1-2,0.5) in (a) and (1.0,0.5) in (b) to match with 
the simulation results; they are of order unity, as expected. The transition anisotropy is 
given by ((J\\/ct ± ) c « [1 + 3n 2 (xp/c A ) 2 ] 1/2 , and (a\\/a±) c ~ 2.5 for (a) and « 2.9 for (b). 
These curves are in good agreement with the simulation results in their respective regions. 
It is evident that there are two saturation regimes. We also observe that the assumption of 
the constancy of ca and xp holds true over a wide range of anisotropy in both the figures, 
although the theoretical curves deviate slightly from those of the simulation results in the 
weak anisotropy side. 

Other comparisons are shown below for electron-positron plasmas. Using Eq. (fTj). we can 
indirectly estimate r from the simulation results: 

r = (2^)-\B/B*)(S/ll)(h ot /h)-%, (12) 

where S is the area of the simulation box and /tot is the total current in one direction along 
the z-axis. Figure 2(a) compares the radius at saturation obtained from this equation with 
those of the model (r A [Eqs. © and JTUD with a( = 0.7] and f P [Eq. ©]). Figure 2(b) 
shows the current per beam at saturation estimated using the radius of Eq. (fT2j) normalized 



by the Alfven current I a = Io&\\, which does not include the isotropization effect. Hence, 
the normalized values are not expected to be unity but to be close to the isotropization 
factor; xa — Ca«C = 0.84 for the Alfven limit, and Xp = 0-5 for the particle limit (shown 
by dotted lines). We observe that the model is consistent with the simulations. Figure 3 
shows the current density at saturation obtained from the simulations for three cases: (a) 
Alfven limit regime (cry = 0.15, x — 0.84), (b) particle limit regime (cry = 3.1, x — 0-5), and 
(c) the transition point (cry = 0.3, x — 0.5). The current density is normalized by Jp, which 
includes the isotropization effect. It is evident that saturation occurred when \J Z \ pa J P for 

(b) and (c), while \ J Z \ <C Jp for (a). It should be noted that in some regions of (b) and 

(c) , the current density exceeds Jp. This occurs because the current-carrying beams are 
pinched. In any case, the typical magnetic field can still be estimated from Eq. ([TJ because 
it is mainly determined by the total current in a beam, which approximates to the Alfven 
current at saturation irrespective of whether the beam is pinched or not. 

Discussion — Even in the presence of a background magnetic field, the proposed model 
is applicable if the time taken for saturation is shorter than the deflection time due to the 
background magnetic field. Otherwise, the problem of the magnetized Weibel instability or 
the whistler instability for electron-proton plasmas ^| should be considered. 

It is shown that the proposed model is consistent with several conditions obtained pre- 
viously. Medvedev and Loeb [^J proposed that saturation occurs when the effective Larmor 
radius ri{B) becomes comparable to the most unstable wavelength in linear theory, Ao- 
If Ao is replaced by the typical beam radius at saturation, f, we obtain the condition of 
^(-Bmax) ~ f, which is qualitatively equivalent to condition (JHJ). Califano et al. lj| also 
found that saturation occurs when tl{B) ~ Iq. This result is in agreement with the result 
of the particle limit in Eq. (JHJ). 

In electron-proton plasmas, if the initial velocity distribution of the protons is the same as 
that of electrons, the upper limits for the proton currents will be given by I' A = (m p /m e )lA 
and I'p = Ip with fj, — 1. Thus, for the particle limit, both fp and -B max would be (m p /m e ) 1//2 
times larger than those of the electron currents. 

In conclusion, the magnetic field generated by the Weibel instability saturates when the 
currents in the beams evolve into the Alfven current; there are two saturation regimes: the 
Alfven limit and the particle limit. The beam model proposed in this letter provides a good 
estimate of the magnetic field strength at saturation. This model will also be useful to 
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consider the evolution of magnetic fields even after saturation. 
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FIG. 1: Maximum magnetic field strength B max for a± = 0.1 as a function of the initial anisotropy, 
cr \\/ cr ± ~ 1 ; ( a ) electron-positron plasma, and (b) electron-proton plasma. Dots denote the results 
of numerical simulations. The dashed and solid curves represent the theoretical results for the 
Alfven limit and particle limit, respectively [see Eq. (|11|)]. This shows good agreement between 
the predictions of the proposed theory and the simulation results in each saturation regime. 
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FIG. 2: (a) Typical beam radius at saturation for the electron-positron plasma as a function 
of the initial anisotropy with a± = 0.1. Dots represent the radii of Eq. ()12|) evaluated using the 
simulation results. The dashed curve represents va [see Eq. (jBJ)]; solid curve 1.2 fp [see Eq. (jHJ)]. 
(b) Current per beam at saturation obtained from the simulation results normalized by the Alfven 
current (I a = Io&\\) ■ The expected values are represented by the dotted lines (see text). The 
results of the theoretical model are consistent with those of the numerical simulations. 
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FIG. 3: Contour plots of current density at saturation, J z , normalized by Jp, on the x — y plane 
for cr_|_ = 0.1: (a) Alfven limit regime (cry = 0.15), (b) particle limit regime (cry = 3.1), and (c) 
the transition anisotropy (cni = 0.3). The spatial unit is the electron skin depth Iq. The dotted 
curves indicate the levels of J z /Jp = ±0.1 in (a), and those of J z /Jp = ±1 in (b) and (c). The 
solid curves indicate the levels of J z = in all figures. We observe that saturation occurs when 
\Jz\ < Jp hi ( a ) an d | J z I ~ Jp in (b) and (c), as predicted by the proposed theory. This figure also 
confirms that the beams are approximately uniform in radius and current for all the cases, which 
are assumed in the theoretical model. 
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